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Abstract 

Biochemical reaction systems in mesoscopic volume, under sustained environmental chem- 
ical gradient(s), can have multiple stochastic attractors. Two distinct mechanisms ai^e known 
for their origins: (a) Stochastic single-molecule events, such as gene expression, with slow 
gene on-off dynamics; and (b) nonlinear networks with feedbacks. These two mechanisms 
yield different volume dependence for the sojourn time of an attractor. As in the classic Ar- 
rhenius theory for temperature dependent transition rates, a landscape perspective provides a 
natural framework for the system's behavior. However, due to the nonequilibrium nature of 
the open chemical systems, the landscape, and the attractors it represents, are all themselves 
emergent properties of complex, mesoscopic dynamics. In terms of the landscape, we show 
a generalization of Kramers' approach is possible to provide a rate theory. The emergence 
of attractors is a form of self-organization in the mesoscopic system; stochastic attractors in 
biochemical systems such as gene regulation and cellular signaling are naturally inheritable 
via cell division. Delbriick-Gillespie's mesoscopic reaction system theory, therefore, provides 
a biochemical basis for spontaneous isogenetic switching and canalization. 

1 Introduction 

Epigenetic inheritance at the cellular level preserves certain phenotypes through cell divisions 
r i[3Tl l. Since the definition of "epigenetic" means the inheritability is not due to genes, i.e.. 
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DNA sequences, the epigenetic phenomenon must be a biochemical process r i[39ll53l l. DNA 
and histone modifications, gene regulations by transcriptional factors, signaling networks, and 
metabolic pathways are all parts of cellular biochemistry. Current research chiefly focuses on 
the "code" of epigenetic inheritance in terms of DNA methylation and/or histone acetylation 
[15611291 [361]. 

While the gene expressions and their regulations are a central component of epigenetic 
processes, it is less certain what the roles of cellular signaling networks, or metabolic net- 
works, ai^e. Even more importantly, ai^e different gene expressions the cause of the epigenetic 
phenomenon, or consequences of the dynamics of a larger intracellular biochemical network 
as a whole? In the present paper, we advance a theory for biochemical reaction systems in a 
mesoscopic volume r B2l l. Taking a broader perspective that is rooted in stochastic, nonlinear 
dynamical systems r B3l l. we illustrate that it is likely that the code for epigenetic inheri- 
tance is distributive [ 1261 1. Using two classes of biochemical networks (see Fig. [T|l as illus- 
trations, we investigate mesoscopic, nonlinear biochemical reaction networks with multiple 
stochastic attractors [ ||60l l. More importantly, we shall show that two very different types of 
mesoscopic bistabilities exist: the stochastic bistability which has no macroscopic counterpart 
[ |[30ll50ll27l[49ll63l l6ll. and the nonlinear bistability which exhibits deterministic bistability 
when the system's volume tends to macroscopic scale [ |l58l[T9l[35ll5Tl l. 

There is a growing awareness of mesoscopic biochemical bistability [ ||37l [1411551 [T2l l9ll. 
Multistability in a mesoscopic chemical reaction system can be mathematically represented 
by the Delbriick-Gillespie processes (DGP) whose probability distribution function follows 
chemical master equation [ |[T3ll 1 and whose exact stochastic trajectories can be obtained by the 
widely known Gillespie algorithm [|[3l[25l]. This approach has found numerous applications 
in recent studies. For an introduction to this chemical reaction system theory, see [ H3l[42l[35l 
[48l[62l[5l]. A review by H6l is particularly accessible. 

In principle, a biochemical reaction system in a small volume of the size of a cell, on a 
very large time scale (years, hundreds years), will have a steady state probability distribution 
which reflects the continuous jumping among the multiple stochastic attractors: regions with 
high local probability [ |[20l[T9l[66l l. In terms of this stationary distribution, the dynamics of a 
mesoscopic biochemical system can be cogently visualized and even quantified by a landscape 
representation [|[T8l[66l[6n[28[[2ni. 

It will be illustrated in this paper that multistability is not only stable against noise, i.e., ro- 
bust, but also could be readily inherited during the process of cell volume change and division. 
Within a cell cycle, it is not necessary to directly control the concentration of a biochemical 
species: As soon as cellular concentrations are deviated from their locally most stable values, 
the dynamics will cause them to spontaneously relax back as long as the deviations are within 
the limit of the basin of the attraction (stable state). Indeed, this restoring phenomenon should 
be regarded as a form of "self-organization". Such state naturally has stability and robustness. 
On the larger landscape scale, the system is "digital". This kind of "inheritable code" is not 
only distributive but also dynamic [ |[26l l. in contrast to Watson-Crick basepairng. 

Some mathematical analysis of the DGP is at the foundation of our current landscape the- 
ory [ B2I [43I I. In Sec. [Hwe offer some discussions. Intuitive and appealing as it is, the 
justification of using the stationary probability distribution as the landscape involves subtle 
mathematical ideas which deserves further investigations [ |[20[[2T| 1. 
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Figure 1 : A biochemical signaling network of phosphorylation-dephosphorylation cycle (PdPC) 
and gene expression regulatory network for a self-regulating gene are essentially isomorphic, (a) 
An enzyme E can be phosphorylated to become E*, catalyzed by a kinase K. Dephosphorylation 
of E* is catalyzed by phosphatase P. The activity of K, however, is modified through binding x 
copies of E*: K^^ denotes K{E*)^ complex, (b) The expression of a transcription factor (TF) is 
regulated by the very TF, in monomer or dimer form (x = 1, 2). The TF is degradated with rate k. 
gi and go are the respective biosynthesis rates when the TF is on and off the gene. If gi < go, the 
TF is a repressor; and if gi > go the TF is an activator. [Redrawn based on Fig. 14 of [|23l .1 

2 Biosynthesis of Self-regulating Repressor with Slow 
On-and-off Gene Fluctuations 

We first consider a simple model for the biosynthesis and degradation of a repressor with 
stochastic gene expression, which is regulated by the repressor. The canonical kinetic scheme 
for this model, neglecting the intermediate stage of mRNA, is in Fig. \Vp. For simplicity, 
we assume that the repressor-gene binding rate h{n) = hoU which involves a monomer, and 
go > gi- See |[30l . |[50l |[27l . and ||60l for extensive studies of this model, and related systems 
with self-activation dimer (/i(n) oc n(n — 1), gi > go). We choose this model to demonstrate 
bistability due to slow, nonadiabatic fluctuations in the gene state. This is a stochastic effect 
due to single-molecule behavior [lH |63l] which disappears in the macroscopic Law of Mass 
Action [||271|43] (see Methods). 

2.1 Stochastic bimodality and bistability 

It is generally accepted that, in biochemistry, noise-induced bistability means a small ki- 
netic system has bimodal distribution while its macroscopic counterpart has only uni-stability 
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Figure 2: Steady state probability distributions for the repressor molecule numbers in mesoscopic 
systems with different volume V but same DNA concentration Xt = I. (A) One gene with V = 1. 
The red curve is the repressor distribution. The two thin curves are p'^'^{n\m) with m = 1,0. Model 
parameters are taken from |[27l . (B) Four genes with V = A. The thin curves are for p'''^(n|m) with 
m = 4, 3, 2. (C) Ten genes with V = 10. (D) Distributions for the concentration of the repressor 
protein as functions of increasing V approaching to the macroscopic prediction with steady state 
at y* = 27.8 (see Methods). (E) The mean exit time for the ail-off, low-transcription state, defined 
when possible as the well at n = Mgi/k = {giXt/k)V, decreases with V . 

r i[27ll49l l6ll. Usually, the meaning of "macroscopic counterpait" is defined as the same chem- 
ical or biochemical reaction systems at same concentrations. For a large volume, the concen- 
trations are deterministic variables; but in a mesoscopic volume, the copy numbers fluctuate. 
This "con^espondence principle" is consistent with the experimental practices: In the past most 
biochemical experiments on gene expression were canied out with DNA measured in concen- 
trations ri[59ll. 

With increasing volume of the a reaction system {V) and the copy numbers DNA (M), 
stability of one of the two states decreases while the other increases. For a sufficiently lai^ge V, 
the bistability disappears all together r il27ll47l l. We call such bistability stochastic bistability. 
Fig. |2] shows how the stationary probability distributions change with the increasing V and 
DNA copy number M while keeping its concentration xt = M /V = 1. 

In agreement with the visual landscapes in Figs. |2l\-D, also shows how the mean 
sojourn time for the ail-off, low transcription state (m = 0, n < 80 in Eq. H]) decreases as a 
function of This is in sharp contrast to the nonlinear bistability (Fig. O we shall discuss 
next. 

'The landscape representation in the nonadiabatic analysis is valid quantitatively: The transition rates between 
any two states a and b, kab and kba, satisfy kab/kta = Pb' /p^^ where p** is the probability of state x. Let's assume 
the shape of each peak is approximately Gaussian. Then the logarithms of the peak values are ln(p^''/A/27rcra) and 
hv{pl^ / \/2TTab), where cr's are the variances of the Gaussian distributions. Thus the "energy difference" between the 
two wells £^(,-£^0 = -ln(fcab/(Th) + ln(fcf,a/o-a), orln(A;a6/fc&a) = ~{Eb-\TLCTb) + {Ea-\TLaa). The right-hand-side 
are the "free energy" difference of the two wells. 
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3 Phosphorylation Dephosphorylation Cycle with Non- 
linear Feedbacks 

A gene regulatory network with dimeric activator also exhibits bistability r il30ll60l[5Tl l. but by 
a different mechanism. To illustrate this, and also to broaden the scope of our discussions, we 
shall consider a cellular phosphorylation dephosphorylation cycle (PdPC) signaling network 
in Fig. [T^. In particular, we shall consider the case of positive feedback with a dimer (x = 2). 
The same analysis can be applied to the gene expression system with x = 2 and gi > go in 
Fig. [TJ). It is easy to show that for x = 0, 1, the cellular signaling network in Fig. [T^, with 
slow fluctuating kinase activity, exhibits stochastic bistability [IHIITI]. The theory we present 
here is general to both types of biochemical networks in terms of nonlinear chemical kinetics. 

The class of networks in Fig. [T^ has been widely implicated, such as in Src family ki- 
nase membrane signaling r ifTTI l. Rab 5 GTPase in endocytic pathway r il65l l. Xenopus oocytes 
regulation for cell fate r ifT5l l. and long-term neural memory r i[37l l9ll. The network has been 
studied in |[T6ll45l[T9l l6l. The detailed kinetic scheme is given in Eq. [Ufor which a DGP is 
uniquely specified (see Methods): 

E + K'< + ATP h E* + + ADP, (1) 

fc_i 

K + 2E*^K\ E* + P^E + P + Pi. 

Note that the stochastic Delbriick-Gillespie approach is not an alternative to the traditional 
enzyme kinetic modeling (Eqs. \T0\ [TT]) . When copy numbers in a chemical reaction system 
are large, a Delbriick-Gillespie process (DGP) automatically yields the deterministic dynamics 
predicted by the traditional Law of Mass Action [||5l|46l]. One of the most important predic- 
tions of the DGP theory is a dynamic landscape for the system, as shown in Fig. [3] as well as 
Figs. |2K-D. Such a landscacpe can only be rigorously defined in a stochastic model; it can be 
computed using a chemical master equation. 

4 Time Scales, Emergent Landscape and Implications 
to Epigenetics 

4.1 Three time scales of cellular dynamics 

Double-well landscapes shown in Figs. |2] and |3] suggest multiple time scales in the biochemical 
dynamics. In fact, there aie three distinct time scales in such systems. Note that while fluc- 
tuations modify the "down-hill" deterministic kinetics, they lead to "up-hill" dynamics which 
is impossible in a macroscopic system. The time for "barrier crossing", however, is extremely 
slow in comparison to the down-hill kinetics. 

Therefore, the fast time scale in the system is the individual biomolecular reactions in Fig. 
[T] For the present work, they are given in terms of the rate parameters fc's in Fig. [T^ and 
/, g, h, k in Fig. \Vp (or equivalently, the a, /?, e, 6 in Eq. |9]) Millisecond are not unreasonable, 
even though individual biochemical reactions inside a cell could be much faster or slower. 
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Figure 3: The landscape for the PdPC kinetic system in Eq. [H 0(a;), according to Eq. [T6l using 
parameters Xt = 1, a = 43, /3 = 10, e = 0.01 and 5 = 0.5. From Eq. [TTl the steady states are 
xl = 0.05 (stable), = 0.632 (stable), and x* = 0.368 (unstable). 

The middle time scale is the relaxation kinetics of a network to its steady states, as illus- 
trated in the Fig. [3]by the downhill dynamics. Note that the very existence of a steady state 
(an attractor), or steady states, is a consequence of "self-organization" of a complex reaction 
network. 

The slow time scale is the transition rates between the two basins of attraction, or "wells". 
Both the middle (deterministic) and slow (stochastic) time scales are emergent properties of 
the biochemical network. Following ||23l we shall denote them molecular signaling time scale 
(MSts), biochemical network time scale (BNts), and cellulai" evolution time scale (CEts), re- 
spectively. 

Again taking the signaling system in Fig. [T^ (also Eq. [Hi as an example. For biochemically 
realistic situations, e ^ a. We thus have the landscape given in Eq. [16] simplified into 



The parameters a, /3, 5 in Eq. [2]define the MSts. We can now use the model to investigate the 
role of a, /3 and 5 on the BNts and CEts. The BNts is given by the ri and r2 in Eq. [121 and the 
CEts is given by the Ti^2 and T2-^i in Eqs. [TSlHOl 

The BNts changes with total concentrations of the regulators E inside the system, as well 
as the concentrations of other factors. Fig. [4] shows how in the simple model the time-scale 
decrease, i.e., rate increases, with the total concentration of E, xt- Within less than one order 
of magnitude change of xt, from 1 to 6//M, the relaxation rate in the state 2, i.e., the well 
on the right in Fig. [3l increase by a factor of 100. Eq. [12] confirms that there is a square 
dependence of the relaxation rate to the concentration. 



(/>(x) = Xt ln(xj — x) — X In 



(ax^ + 5){xt — x) 




(2) 
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Figure 4: Relaxation kinetics in the state 2, the well on the right in Fig. [3] corresponding to state 
as a function of the total concentration of E, Xt, according to Eq. [121 The parameter used are 

a = 4:.3fj,sec~^fiM~^, (3 = lOmsec"^, e = lmsec~^ fiM'"^ and 6 = 0.5msec~^. 

The CEts is extremely sensitive to the number of molecules in the biochemical system. Fig. 
|5] shows that with the MSts on the order of milli- and microsecond, and with concentration of 
E on the order of micromolar , the transition times between the two states in Fig. |3]can be as 
long as thirty thousand years! Hence, the stability of the emergent attractors could be extremely 
robust against spontaneous concentration fluctuations (i.e., intrinsic noise) in the system. We 
also notice that both transition rates decrease with the V. This will be explained in Sec. 14.31 
below. 

More interesting biologically, we note that the range of 700-1000 copy number of E cor- 
responds to a time range of 10 hours to 30 years. In yeast, ll24l have reported that the copy 
numbers for most of the transcription factors are centered arround 2^*^ = 1024 per cell. 

4.2 Epigenetic inheritance and canalization on the CEts 

Let us consider two replica of a mesoscopic biochemical reaction system in a laboratory, for 
example one of those in Fig. [T^ which do not involve gene expression. If the two systems have 
same total E but different initial values for E* , one near zero and one near the total E, then 
these two systems settle to the two different attractors. In the time scale much shorter than 
the evolutionary transitions, the numbers of E* fluctuate around the and 722, respectively, 
or equivalently around the c\ and C2 if the volume of systems do not change. (If the volumes 
are changing, then the fluctuations are ai^ound the and C2, but not and 722!) However, 
at the time scale greater than the cellular evolution, there will be transitions between the two 
attractors. This is shown in Fig. [6]\. The probability distribution for the concentration of E* 
is shown in Fig. |6]C. 

Fig. [6^ shows the identical reaction system, except its volume and total E are twice as 
large (keeping the concentration invariant). What we observe from the Fig. |6ti; is that the 
"concentration" distribution for E* in the two cases have essentially the same locations for the 
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Figure 5: The switching times between two stable states, Ti^2 (dashed line) and T2^i (solid 
line) increase with the volume V in nonlinear bistability (according to Eqs. [TSl) . Xt = 0.1 /iM, 
a = 4.3/isec~^, (3 = 10msec~\ e = lmsec~^ and S = O.Smsec"^. A cell has a volume at about 
10-20 femtoliter. With concentration xt = 0.1 /iM, there are 600-1200 number of total molecules. 
This gives the switching time around 3 x 10^ to 10^^ sec. That is 10 hours to thirty thousand years! 
Therefore, with rather rapid individual biochemical reaction rates a, (3, e and 5, the emergent 
epigenetic states can be extremely stable. 

peaks and trough. This implies that if the size of the biochemical reaction system increases in 
the time scale sufficiently short, then the identities of the attractors can be preserved. 

The stable state of the system is not only stable against intrinsic noise, but also could 
be readily transfeiTcd to the two daughter cells. During the cell cycle, the concentrations of 
biochemical substances might become approximately one half of the original value, in the 
extreme cases due to cell volume increase, and the system deviates from its corresponding 
stable state. The kinetic law of these two daughter systems is just the same as their mother 
cell except they have not relaxed to any of the stable states. As we have shown previously 
that the relaxation scale is not very large, so they will spontaneously relax back to the same 
corresponding stable state as long as not leaving an basin. Many previous work of epigenetics 
always searched for a stable chemical substance like DNA, which could be self-propagated 
and inherited to the daughter cells, while here we give another alternative possibility that the 
code of epigenetic is at the dynamic level of the whole biochemical network. 

Further, there ai^e clear upper and lower bounds for the rate of volume increase: it can 
not be too large such that the instantaneous changing concentration is outside the basin of an 
attractor; it can not be too slow such that it is on the cellular evolution time scale. Surely, the 
stability of epigenetic code is weaker than that for a stable chemical substance and it is more 
flexible facing the influence of the fluctuating environment, but it is suffient for a normal cell, 
a chemical system, to survive and inherit even in a fluctuating environment. 

Fig. |6]D shows precisely two of such simulations. Consider the volume, and the total E, 
double within the time of 5 units. This is a duration much shorter than the cellular evolution 
time. The top and bottom traces are the number of E* from two simulations. At the end of 
the doubling, if each system is divided into two, both "daughter systems" will also inherit the 
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Figure 6: Nonequilibrium steady state fluctuation of the number of phosphorylated E* as functions 
of time, with V = 100 in (A) and V = 200 in (B). The steady state distributions of the number 
of E* for the V = 100, ut = 100 and V = 200, ut = 200 are shown in (C). The parameters used 
are a = 43, /3 = 10, e = 0.01 and 6 = 0.5. (D) The volume V and rit are increased by a factor of 
2 from V = 100, rit = 100 to V = 200, rit = 200 within the time to 5. One sees that the two 
attractors are well separated in the volume doubling process. 

biochemical state of the "mother system". The biochemistry of a mesoscopic reaction system 
is inheritable! The epigenetic stability, i.e., canalization r il36l l could be directly related to the 
CEts. 

4.3 Mesoscopic stochastic dynamics on a landscape 

The foregoing discussion clearly illustrates the power of the "landscape" perspective in visu- 
alizing and representing the global dynamics in bistable systems. We see that for stochastic 
bistability, at least one of the "barrier heights" decrease with volume V (Fig. while for 
nonlinear bistability, the both barrier heights increase with the volume (Figs. |5]and|6ll. 

What determines the entire landscape? Since it is defined on the space for all possible 
concentrations and/or copy numbers of all the moleculai^ species in the reaction system, it- 
self can not be determined by the concentrations and copy numbers. Rather, it is determined 
by the all possible molecules involved and their interaction/reaction rate constants. In other 
words, biochemical reaction networks. Since the molecular interaction/reaction rate constants 
are properties of molecular structures, which in turn is determined by the primary sequences 
in the case of proteins, we conclude that the landscape, conceptually, is encoded in the DNA 
sequence, together with the extracellular environment including the volume V, but is inde- 
pendent of the expression patterns of transcription factors. They are the consequences of a 
biochemical reaction system's dynamics [[50]]. 

There are many similaiities between the energy landscape for a protein in equilibrium 
[HH] and the landscape for an open mesoscopic chemical system in a nonequilibrium steady 
state r i[66ll4ni23l l. We, however, want to emphasize a key difference: Recall that an energy 
landscape exists a priori for a dynamical protein tl8i . The open-chemical systems are fun- 
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damentally different in this respect r POl 1411 |23ll. Specifically, for any finite volume V, a 
mesoscopic reaction system has a stationary probability distribution for the number of copies 
of all its biochemical species, pv'(n), where n = (nx, ray , n^, • • • ) are the copy numbers of 
the molecules X, Y, Z etc. It can be shown that such a distribution can be expressed as 



py(n) = exp 



-y^(x) + ^,(x) + ^ 



(3) 



where x = n/V. Furthermore, it can be shown that the function <j){^) is a meaningful land- 
scape for the complex dynamics of the nonlinear biochemical system. That is, the dynamics 
always goes "down the hill" of 0(x) [ 112811211 1. though usually not by the steepest descent path. 
As we have seen, the landscape provides an very useful organizational device for thinking 
about cellular biochemical dynamics at widely different time scales, ranging from individual 
signaling reactions to cellular phenotype switching. 

Eq. |3] shows that the stationary distribution py(x) changes with V. When the V becomes 
macroscopic volume, the probabilities are concentrated only at the global minima of the func- 
tion (/>(x). However, with increasing V, — {1/V) Inpv (n) approach to the function (/)(x) which 
is defined on the entire space of x. It is an important insight that such a landscape exist and it 
is independent of the system's volume, as long as the volume is reasonably large r il28ll2Tl l. 

The nonlinear bistability, therefore, is the mesoscopic manifestation of a double well in the 
(/)(x). Its macroscopic counterpart has two stable steady states. However, stochastic bistability 
is something very different: The bimodal distribution only exists when the V is sufficiently 
small; when (pi and cl)2 in Eq. [3] contribute to the py(n). 0(x) has only a single well. With 
increasing volume: the barrier in the nonlinear bistability increases, while that of stochastic 
bistability decreases. 

The emergent landscape of cellular interaction network dynamics and the landscape for 
protein dynamics are fundamentally different. The lack of detailed balance due to the open 
chemical nature of the former gives rise to the cycle flux underneath the landscape r i[6Tll23l l. 
The cycle flux makes the landscape non-local. When such a flux is sufficiently strong (i.e., 
mathematically characterized by the emerging of complex eigenvalue and eigenvectors), a syn- 
chronized dynamics arises r il44ll22l l. The emergence of synchronized dynamics requires an 
entirely new kind of phenomenology. In the macroscopic classical world, this is the birth of 
oscillatory behavior and wave phenomena that have ruled classical engineering for a century. 



5 Discussion 

5.1 Adaptive landscape 

While the concept of energy landscape becoming a very useful term in protein dynamics, 
the concept of adaptive landscape in evolutionary dynamics is still highly controversial [HI]. 
We believe one of the main reasons for this situation is that the landscape in the latter, as 
the landscape in the present work, is an emergent entity, which is not given a priori. There 
are two issues related to this important distinction: (1) The mathematical existence of such a 
landscape in a general, nonlinear stochastic dynamics which does not have detailed balance; 
and (2) How is such a landscape, even exists, related to the dynamics, both deterministic and 
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stochastic. The most nontrivial issue here is the logical relationship between the landscape 
and the dynamics: It is rather clear that in systems with detailed balance, the landscape exists 
a priori and the dynamics is a consequence of the landscape. However, for system without 
detailed balance, the dynamics, as defined by the reaction networks and all the individual 
rate constants, define the overall dynamics as well as define the landscape. Dynamics and 
the landscape have correlations but no causality. Hence, logically, it is not con^ect to view 
the dynamics as a consequence of the landscape. Nevertheless, the landscape is still a very 
useful device to understand and characterize the overall dynamics. The concept of emergent 
landscape, thus, should only be understood in this "historically retrospective" sense. 

5.2 Genocentric epigenetic inheritable memory and a possible al- 
ternative 

Currently, methylation of DNA is considered to be the leading candidate for epigenetic inher- 
itance. Even though this mechanism is not based on Watson-Crick basepairing, its function 
is still intimately dependent on the discovery of 1950s. Namely, the "code" is still a part of 
an extended, modified DNA structure. While certainly methylation is a part of the "whole 
picture" of epigenetic regulation, such a view might be too genocentric and too limited. It is 
not unlikely that the code of epigenetic inheritance is dynamics rather than static, distributive 
rather than localized: It is an emergent property of the whole system rather than depend on a 
very few number of substances or regulatory mechanisms. The emergent landscape of a meso- 
scopic chemical reaction system share certain features with the content addressable memory 
proposed by J.J. Hopfield years ago, but eliminated the technical need for detailed balance 

We would like to point out the fundamental difference of our proposal is a non-genome, 
pure biochemical based epigenetic mechanism. While DNA methylation is a part of this pic- 
ture, our proposed mechanism moves the focus away from DNA basepair recognition and 
memory, and shift it to biochemical networks. The memory in our model, in principle, can be 
independent of DNA. We understand that immunological diversity and memory has now been 
shown to be DNA based r il54l l. Still, the digital "quantal" nature of immunity at individual 
cell level has been noted r il52l l. Also, in the field of neuroplasticity, cellular mechanism for 
memory has been proposed to be very similar to our PdPC model r i[37l l9ll. A biochemical 
based memory is too natural to be completely neglected by evolution. 

At single cell level, bimodality has been regularly observed in diverse cell types. Recently, 
for example, ifTOl demonstrated in their Fig. IC bimodal expression level of lactose permease 
in E. coli. corresponding to the two states of a bacteria cell, uninduced or induced by lactose 
analog TMG (methyl-b-D-thiogalactoside). STJ in their Figure 4 reported bimodal Raman 
spectra as an indicator for DNA fragmentation in apoptosis of DAOY cell line (human brain 
tumor meduUoblatoma), and in Figures 6 and 7, |[64l showed bimodal FITC-Annexin V protein 
in apoptosis of U20S cell line (human osteosarcoma). 

5.3 What is a pathway, and what are cross-talks? 

The concept of biochemical regulatory pathway is widely accepted in the molecular cellular bi- 
ology. In general apphcations, a pathway provides a sequential events of activation/deactivation 
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in terms of the regulatory/signaling proteins. However, from our biochemical reaction system 
perspective, the state of a cell, as a stochastic attractor, is defined by the states of all the regula- 
tory/signaling proteins [(221]. Hence, from a more rigorous theoretical standpoint, one needs 
to know, when a protein Y is activated, what are the states of its upper-stream and down- 
stream proteins, X and Z. In this sense, the "pathway view" of cellular signaling is similar- to 
the "structural pathway" view of protein conformational change. While it is useful, it can be 
mis-leading r ifTTl l. 

Almost universally in the discussion of cellular signaling, the concept of "certain pathway 
leading to certain response" is an established language. But knowing a great deal of "cross- 
talks" in signaling pathways r i[34l [38l l. this linear thinking is incorrect. In fact, a response 
is really a changing state of a cell, be it proliferation versus growth arrest, apoptosis versus 
senescence. So they should not be associated with only one pathway or another. Rather, they 
are different states of entire integrated network. It is true, by "activating certain pathway" 
while keeping other the same, one might promote a particular state, but that does not imply the 
response is only associated with that pathway. 

Now if we take the view these responses are different "states" of a cell, then many path- 
ways are involved. And the most important things about biological complexity is how many 
these stable states are "available" in the system r il50l l. More states also means more possible 
transitions between them (whether a transition actually occurs is an issue of time scales), and 
thus the complexity is associated with multi-stability. This view is consistent with the defini- 
tion of cellular complexity by possible number of responses to combinatorial stimuli r il38l l. 
If there is only one state, then the system always goes back to somewhere it starts, then no 
complexity. 

The key idea here is not to associate a particular cellular response to a particular pathway; 
but rather activating particular pathway promotes certain response which is defined as state 
change. 

5.4 Living matter: Mesoscopic open chemical system as a bio- 
chemical machine 

ll33l have discussed the possible new phenomena at mesoscopic scale which they called the 
middle way. Others has asked "what is and how is living matter different from the three (or 
five) known states of matters from physics [ID]. The study of mesosopic chemical systems 
offers some insights: 

1) If a chemical system is too large (for example, grinding all the cells into a single tube 
without the small volume of a cell), then the chemistry is different. It completely loses the pos- 
sibility of multi-stability at that CEts [ [llQiH . which is a defining feature of complex dynamics 
[ ||47l ]. The size of biological cells might indeed be a consequnce of living matters possessing 
mesoscopic complexity. 

2) This idea is generally understood but its implication has not been widely appreciated: 
A living matter has to be in continuous exchange of materials, with chemical gradient, with 
its environment; The driven nature of a living matter is completely different from the classical 
way of thinking a "matter" which is being in isolation. To define a living matter, one can not 
completely separate the system from its environment: This is the origin of systems view of 
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holistic biology. 

3) The classical thesis of "irreversibility" of Boltzmann is to understand the spontaneous 
processes leading to equilibrium. This is a different problem as that suggested by Schrodinger 
in his "what is life". The former is to understand the macroscopic in^eversibility in a system 
with Newtonian mechanics, while the latter is a completely different problem; Life is an open 
system: To a first-order approximation, it is an isothermal system sustained by an active input 
and output of chemicals with Gibbs free energy difference. The irreversibility of such a system 
is self-evident. Hence, essential thesis in nonequilibrium physics of living matter is not that of 
Boltzmann, but more quantitative understanding of how such chemically driven systems give 
rise to complex "living" behavior such as self-organization and inheritability. 

4) The significance of the bistability arising from the Delbriick-Gillespie is not that it has 
two different "phenotypical" states for relatively low and relatively high levels of "inducer", 
borrowing the language from Lac operon r ifTOl l. but that both states co-exist for an range of 
intermediate level of inducer! This is reflected by the distribution in this intermediate range 
being bimodal. Such a realization was emphasized in protein folding by [32|. This realization 
can have important implications: A pre-cancerous state might already exist "on the other side 
of the mountain", which is encoded in our genome ri[2ll. 



6 Methods 

6.1 Self-regulating gene and stochastic bistability 

We consider the coupled birth-death process for the gene regulatory network in Fig. [TJ5 with 
X = 1. We note that in a macroscopic biochemical experiment the rate of protein synthesis per 
DNA is (i = 0, 1). Furthermore, we denote the concentration of the DNA with repressor 
X = m/V and the concentration of the repressor y = n/V . Then 

dp{m, n) 



dt 

= g{m)p{m,n — 1) + {n + l)kp{m,n + 1) 
h ft 

+ -^{M -m-l)p{m-l,n) (4) 

5((m) + nk + ^^^(M — m) + mf^ p{rn, n) 

+ (m + l)/p(?Ti + 1, n), 

where g{m) = go{M — m) + gim, M = xtV is the copy number of DNA. This model, in the 
limit of y — )• 00, yields the Mass Action kinetic equation 

— = hoy[xt- X) - fx, — = go{xt- X) + gix - ky. (5) 

Eq. [5] has two roots but one in the interval (0,Xf). For parameters g^ = 80,^1 = l,k = 
1, ho = 0.007, / = 0.1, k = 1.0, xt = 1, we have the steady state x* = 0.66 and y* = 27.8. 

Using nonadiabatic approximation, we have p{n\m) being a Poisson distribution with 
mean g{m)/k. Then the 2-d model is reduced to 1-d with birth rates ho{M — m)g{m)/{kV) 
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and death rate mf. The stationary distribution for the 1-d model can be analytically studied 
following ||6l. 

Mathematically, one can understand the problem as an eigenvalue perturbation: The Markov 
operator involved has the block structure [[47]]: 

Lo-h f \ 
h Li-fj' 

in which h and f are small and treated as a perturbation. The unperturbed operator has a 
degenerated eigenvalue zero, with eigenvectors on the left (1,1) and (1,-1), and on the right 
(POiPi)^ and (po, — Pi)^- Po and pi are Poisson distributions with mean no = go/k and 
ni = gi/k, which are the stationary distributions for the nonperturbed problem. Now for the 
perturbed system, it is easy to verify that zero is still an eigenvalue Aq = 0; however, there is 
also a nonzero, smallest eigenvalue 



Furthermore, the approximated right eigenvectors associated with Aq and Ai aie precisely 



, _ , . Po, , _ , . Pi ^^'^ (Po,-Pi) 



(8) 



A more accurate approximations for the two eigenvectors can be obtained, if needed, by car- 
rying through the first-order perturbation calculations [ ||5T| 1. The eigenvectors for the Ai has a 
nodal decomposition which yields two connected domains with positive and negative valuesjl 
This provides a rigorous mathematical definition for the two states of the system. Since all the 
other eigenvalues Aj » Ai (i > 2), the dynamics within each of the two domains are on a 
different time scale, and fast equilibrated. The remaining slow dynamics corresponds precisely 
to a two state system [ ||58l l with transition rates h and /. The nonadiabaticity plays a decisive 
role in this problemJl 



6.2 Chemical master equation and landscape 0(a:) 

The PdPC with feedback in Eq. [T] is the same as that in Fig. [T^ with /cq = 0. With rapid 
binding K + 2E* ^ K"^ and assuming » [E*]'^, we have the following autocatalytic, 
nonlinear- chemical reaction system: 

E + 2E*^'iE\ E*^E, (9) 

e 5 

^The mathematical theory for the bistability proceeds with the following argument: If one treats the Ai as a small 
parameter, then for system with Ai = 0, the stochastic dynamics can be reduced to two independent subsystems, each 
with a unique stationary distribution. Therefore, with the error on the order of Ai, both eigenvectors for the Aq = 
and Ai 7^ are linear combinations of the two stationary distributions defined on the subspaces. This explains the 
origin of the two-state behavior on the slow time scale (CEts) ISSl . 

^The problem appears very similar to the quantum mechanical eigenvalue perturbation with degeneracy. However, 
it is worth pointing out that in quantum mechanics one computes the eigenvalues which serves as the "energy land- 
scape" in Heitler-London theory; here we compute the eigenvectors as the "energy landscape". This distinction has 
been noted by Ii47l . 
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in which a = ^[K][ATP], (3 = h[P], e 



[K][ADP], and 5 = k-3[P][Pi], 



where [K] and [P] are the concentrations of the kinase and the phosphatase. One can see that 
this kinetic model is intimately related to the Schlogl model which has been the prototype of 
nonlinear- chemical bistability ri[5ll58ll. 

Let X be the concentration of E* , then the macroscopic kinetic equation for Q in terms of 
the Law of Mass Action is 



dx 
Itt 



ax'^{xt — x) — (3x — ex^ + 6{xt — x), 



(10) 



where xt is the total concentration of the E and E* , assumed to be constant in the system. The 
dynamic has three positive fixed points. For most biochemical applications e <^ a and (5^/3. 
Hence, one can approximately have the three steady states 



Sxt 



Xt + \Jx1- API a 



and 



Xt - \/x1 - 4/3 /a 



(11) 



Hence, when aXf > 4/3, there is bistability. x^ and X2 are stable steady states, and is 
unstable. The basins of attraction for x\ and X2 are [0, X3) and (x3,oo), respectively. The 



corresponding linear relaxation rate for the steady state x*, {i 
2axtx* + (3. That is, 

ri = /3 and r2 = 0x2(2x2 — xt). 
The chemical master equation for the system in Eq. |9]is 

dp{n) 



1,2), is r, = 3a{x*y 



(12) 



dt 
a 



( 



+ 



(n — l)(n — 2) + (5^ (nt — n + l)p{n 
P + ^n{n - 1)) (n + l)p(n + 1) 
1) + {nt - n) 



Vy2 



n n 



+ (/3 + :^(n-l)(n-2))n p{n), 



(13) 



where rij is the total number of E and E* molecules, and V is the volume of the mesoscopic 
system. The steady state probability distribution for the number of E* is 



k=2 



[ak{k - I) + 5V\nt - k) 
[PV^ + ek{k - I)] {k + l)' 



(14) 



where C is a normalization factor. If V and nt are large, but their ratio is hold constant, 
then one can develop an approximated formula for the probability density function for the 
concentration x = n/V: 

r{x) oc e-^'^(^), (15) 
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where the (hix) 



(x) = Xt ln(xt — x) — xln 



(ax^ + 6){xt — x) 
(/3 + ex'^)x 



(16) 



2\ — arctan ( ^ / —x ] + 2\l — arctan ( ^ / —x 
V a 

One can check the extrema of (/){x) by setting its derivative being zero: 



d(j){x) 
dx 



In 



(ax^ + 5){xt — x) 
+ €x'^)x 



0. 



(17) 



We see that the extrema of are precisely the roots of Eq. [TOl x^, x\ and x\: The ratio in 
the logarithm being 1 con^esponds to the right-hand-side of Eq. [lO]being 0. 



6.3 The switching time 

The mean time for switching from x\ attractor to x*2 attractor can be analytically computed for 
1-d model. For discrete case, this formulae has been widely used in the various lattice hopping 
models with birth-death processes ll57l : 



J. ' 1 

r?i=n]|+l 



n=0 



+ E (IS) 

n=nj+l m=n+l 



where n* = [x*^V], and 



Un =[-^{n- l)(n -2) + 5){nt-n + 1), 

Wn={li + ^{n-l){n-2)) n. 

When V ^ oo,Wm w (y) F where 'u;(x) = /3x + ex^, andp'*'*(n) ~ (C/y)e"^'^(^), we 
have (see Appendix) 

2T^pV{<t>(xi)-4>{xi)) 
Ti^2 ~ (19) 

wix*3)V-^"ix*l)^"ixl) 

and similarly, 

2„gV(<A(4)-'!^(^2)) 

^2^1 ~ (20) 

«^(4)v-</'"(3;2)<?^"(4) 

This result, as expected, is very similar to &amers formula for overdamped barrier crossing 
in the energy landscape cl){x) with "frictional coefficient" being l/iu{x'^). Note that at Xg, 
w{x'^) = ^(xg); hence a symmetric expression for the "fractional coefficient" is 2/ {w{x'^) + 
uixD). 
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A Mean Switching Time for 1-D CME 

The mean time for switching from xl attractor to X2 attractor can be analytically computed. 
For discrete case, according to |[57l . let 

Un = [^{n- l)(n -2) + 5^{nt-n + 1), 
Wn= (/3 + :^(n- l)(n-2)) n, 

we have 

+ E ^'"(^) E — lir-y ^21) 

n=x*V+l m=n+l "'^ ^ ^ 



Wm ~ w (y) y, where w{x) = j3x + ex^, andp'''*(n) « {C/V)e ^"^(v), we have 



tl_>2 ~ 

n=0 



+ E y E nn\e-^'t>{f) 

n=xlV+l m.=n+l^\V)^ 

Ve-^'^(^) / = dx 



n=0 

x^V-l 



n=x*y+l 



Jo Jxi w{x)e-y^(^r^ 

+V / e-^<^(3') / ——^^dxdy. 
Jxi Jy wix)e-y^i-) 



Applying Laplace's method, one has 

n e-^'^^'^^dx 
Jo 



w{x)e-V^i-) ^ v/-y</,"(4)' 
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Then 



l->2 



+v 



+v 



1 

D(y)e-^'^(i') 

e 



x\<y <x_ 



3' 



X3 < y < X2. 



(22) 



_ wiy)e-y't>(v) 



Furthermore, 



hence 



-V<t>iy) 



dy 



V<P'{y) 



dy. 



2^VcP"{xl) 



2^gy(<7i(x|)-9i(xJ)) 



w{xl)yj -<j)" {x\)f {xD 



w{y)(j)'{y) 



dy, 



in which the second term is comparable to f^^ u(y)~w{y) which is the time relaxing from 
the unstable fixed point X3 to X2- It should be neglectable when the noise is present. 
Finally, when y — )• 00, we have the switching time in terms of cl){x): 



-l->2 



and 



2^gV(0(x*)-0(x*)) 

w{xl)^-(P"{xl)(P"{xl)'' 

2^gVi4>(x*)-<t>ixi)) 



w{x*)^/-(l)"{x*)(l)"{xl) 

Comparison with Kramers' theory 

The Kramers theory considers a diffusing particle that obeys the Fokker-Planck equation 

dp{x, t) 



(23) 



(24) 



dt 



lj[-^P'{x)pix,t)] + kBT^p{x,t) 



For energy function U (x) with two energy wells, say well a and well b, he derived |27| the 
rate constant for transition from a to c by crossing barrier c, known as the celebrated Kramers' 
formula: 



k. 



27rr7 



£g(C/a-C/c)/fcsT 



(25) 
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where fia and fic are the curvatures of the energy function U{x) at a, the well, and c, the 
transition state. Similarly one has kh^a- Then their ratio, which relates the populations in the 
two wells at equilibrium: 




where Fa = Ua + ksT In ^/JI^ is the free energy of energy well a. 



Comparison with diffusion approximation to CME 



The result in Eq. |23l as far as we know, is new. In the past, analysis of baiTier crossing in the 
CME, in the limit of large V , have been based on diffusion approximation to the CME. In that 
approach, one derives a one-dimensional Fokker-Planck equation for large V: 

dp(x,t) d , s\ f \ 

= -^(^(^) - w{x))p{x,t} 

Here, it suggests a potential 

•*W = -2r44-^d., (27, 

J u(x)^w{x) 

Hx) = - flog^dx. (28) 
J w{x) 

(j){x) is the coiTcct one while ip{x) could give incorrect result for the Ti^2 and T2_j>i. See 
r il58l l for an extensive discussion. 



which is different from our 



B Transition Time in Nonadiabatic Gene Switching 

We consider the coupled birth-death process with probability distribution [pQ (n) , pi (n)] , where 
and 1 represent the gene states, with and without bound transcription factor, and n represent 
the copy number of the transcription factor which is the gene product. The distribution satisfies 
the chemical master equation 

= goPo{n - 1) - (50 + nk) po{n) 

+ {n + l)kpo{n + l) - h{n)po{n) + fpi{n) 

dpi (n) , N / , N / N 

= 9iPo{n - 1) - [gi + nk) po{n) 

+ (n + l)kpo{n + 1) + h{n)po{n) - fpi{n). 
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We shall consider h{n) = hn{n — l)/2. If the gene switching between and 1 is slow, 
i.e., it is nonadiabatic, then we have the rapid pre-steady states (conditional probability) in the 
state i (= 0, 1) following their respective Poisson distribution 

p{n\i) = ^ (^Y e-3^/\ (29) 



nl \ k 

Then the mean transition rate from state to 1, and 1 to 0, ai"e 



ko^, = J2h{n)p{n\0) = ^, ki^o = J2fPi^\^) = f- (30) 

n=0 n=0 

Perturbation method for eigenvalue problem 

We can also solve the eigenvalue problem using the method of perturbation theory. Note that 
the unperturbed system, i.e., when h = f = 0, has degeneracy. For the perturbed system, 
we still have an eigenvalue 0, corresponding to the stationary distribution [p(n|0),p(n|l)]^. 
The other eigenvector on the left is [1, —1] and on the right is [p(n|0), — p(n|l)]^, with the 
corresponding nonzero eigenvalue: 

Ai = -/-g, (31) 

which is exactly the kQ^i + ki^Q in Eq. [30l Expressing the Ai in terms of the nondimension- 
alized parameters in r il60l[5ll l. we have 

_ Ai _ h / hgl f goV 

k ~ 2k + 2fc3 \gj 



1 \ 2 / \ 2 



91 J \a\ 



0.29k. (32) 



The stochastic separatrix is well defined in this case by the domains of positive and negative 
values of the eigenvector associated with eigenvalue Ai r i[58l l. 



Nonadiabatic reduction to 1-d with M copies of DNA 

The 1-d model has the transition rates: 

/io(go(M - m) +5im) 

m — ^ m + 1 : ; (M — m), 

Vk 

m + 1 — 7- m : /(m + 1). (33) 

If we let X = m/y then we have 

2 I ( n I 2 



(fo - 9\)^ + \ 9\^t - ^g^xt ~ J^j ^ + 9oXt = 0. (34) 
This should be compared with the quadratic equation for the model in O : 

{ki + k^i)x'^ - {kixt - k2- k^2)x - k^2Xt = 0. (35) 
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